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We present the two-loop QCD amplitude for the interaction of two gluons and a CP-even Higgs 
boson in the Minimal Supersymmetric Standard Model. We apply a novel numerical method for 
the evaluation of Feynman diagrams with infrared, ultraviolet and threshold singularities. We 
discuss subtleties in the ultraviolet renormalization of the amplitude with conventional dimensional 
regularization, dimensional reduction, and the four dimensional helicity scheme. Finally, we show 
numerical results for scenarios of supersymmetry breaking with a rather challenging phenomenology 
in which the Ifiggs signal in the MSSM is suppressed in comparison to the Standard Model. 



The loop-mediated interaction of a Higgs boson and 
gluons is the main production mechanism for a Higgs bo- 
son at hadron colliders. Depending on the decay channels 
in which a Standard Model Higgs boson may be discov- 
ered, the measurement of the signal cross-sections can be 
achieved with a precision of about ±10% or better 
Amid the discovery of a Higgs boson, the signal cross- 
section will be an independent precision test of the Stan- 
dard Model and its extensions. 

The gluon-fusion Higgs boson production cross-section 
is sensitive to higher order QCD corrections 0, S 0| • 
the LHC, the signal cross-sections are known to change 
up to a factor of two when perturbative corrections from 
O (q!^) through O {ct^ are included The accuracy 

of these theoretical predictions is about ±10%. The par- 
tonic decay width to gluons also increases by a factor of 
two when including the known higher order QCD cor- 
rections. This quantity is now computed through order 
O (ct^s) with an accuracy better than 1% Theoretical 
uncertainties in the ggh interaction within the Standard 
Model are currently adequately small for a future com- 
parison with LHC data at a 10% precision level. 

Extensions of the Standard Model (SM) postulate di- 
verse mechanisms for breaking the electroweak symme- 
try. A different Higgs boson sector than the one in the 
Standard Model and additional new particles are often 
introduced. The effects of undiscovered particles on the 
gluon-fusion cross-section are rather unconstrained from 
experimental data. Novel colored particles can change 
the Higgs and gluon interaction dramatically. For ex- 
ample, an additional heavy-quark with SM-like Yukawa 
coupling to the Higgs boson would contribute almost as 
much as the top-quark to the ggh amplitude. Recent 
examples in well motivated models were shown in Q. 

The complexity of the two-loop SM computations at 
O (a^) in the full theory and at O (af) in the limit of 
a heavy top-quark is serious. The methods that have 
been employed are powerful and may be employed in 
other models. In especially simple modifications of the 
SM Lagrangian, for example adding a fourth generation 
with heavy leptons and quarks, the existing calculations 
are already sufficient. However, it will be important to 



know the gluon-fusion cross-section through at least or- 
der O (ofj) in all models which aspire to explain LHC 
data. Many viable extensions of the SM contain colored 
particles where, at two-loops, more than one of these 
massive particles appears in diagrams contributing to the 
giuon-giuon-higgs amplitude. Known analytic methods 
for two-loop calculations are restricted to problems with 
a small number of mass parameters. New techniques to 
evaluate multi-loop integrals with diverse mass-scales are 
indeed required. 

In this article, we compute the two-loop QCD ampli- 
tude gg ^ h,H for a CP-even light {h) and heavy (H) 
Higgs boson in the minimal supersymmetric extension of 
the Standard Model (MSSM). We employ a numerical 
method which we have recently developed for multi-loop 
calculations [l^. In our method, Feynman diagrams 
with diverse combinations of massive and massless prop- 
agators are treated on an equal footing. Given the com- 
plexity of the MSSM, the computation of 2 ^ 1 two-loop 
amplitudes in other extensions of the SM may also be 
tractable. 

Partial SM-like contributions from supersymmetric di- 
agrams with only squarks in the loo ps h ave been recently 
computed in the literature [H [H, 111, [ll • The MSSM 
amplitude is also known in the limit of a light Hi ggs b oson 
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with respect to quarks, squarks, and gluinos 

The computation of the two-loop amplitude without 
using an effective theory approach is well motivated in 
SM extensions. Relatively light colored particles (lighter 
than the the top-quark) are not excluded experimentally. 
Heavy Higgs bosons are also predicted in the spectrum 
of new theories. Also, the couplings of the bottom quark 
to Higgs bosons may be significantly enhanced in mod- 
els with more than one Higgs doublet. The MSSM ex- 
hibits all of these features; knowledge of the amplitude for 
the gluon-gluon-Higgs interactions without assumptions 
about the mass hierarchy of Higgs bosons and colored 
particles is therefore especially important. We consider 
the MSSM an archetype for many other models regarding 
its computational challenges. To the best of our knowl- 
edge, we present here the first complete result for a two- 
loop three point Green's function in the MSSM. 
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FIG. 1: Sample of two-loop diagrams containing up to four 
different mass parameters in the propagators. 

The gg ^ h amplitude at O (a^) includes 135 two- 
loop Feynman diagrams. We have generated them us- 
ing QGRAF la,]. We implemented the MSSM Feyn- 
man rules following the method of Ref. [l^l for Majorana 
fermions. Traces of Dirac matrices and the color alge- 
bra are carried out with programs written in standard 
algebraic manipulation packages, including FORM [20| . 
We have checked that our computer programs and Fey- 
nArts [21] generate equivalent integrands for the ampli- 
tude. 

We depict some of the Feynman diagrams in Fig. [1] 
We computed diagrams with only one massive particle in 
the loops both analytically [ll| and numerically @. We 
computed the remaining diagrams numerically applying 
the method of [l3] • We review here only the salient 
features of the numerical method. 

We combine propagators in each loop with a set of 
Feynman parameters and integrate out the loop mo- 
menta. The Feynman parameters are mapped to a unit 
hypercube integration domain. At a second step, we ap- 
ply a sector decomposition algorithm 22, 2^] to factorize 
divergences in the Feynman parameters which are reg- 
ulated by the dimensional parameter D = 4 — 2e. We 
deform the contour of integration away from the poles 
which may occur for physical values of the external mo- 
menta and mass parameters (24| . Finally, we extract the 
factorized singularities in e with independent subtrac- 
tions on the real part of Feynman parameters and per- 
form an expansion around e = [§] . With this procedure 
the amplitude is written in terms of 0(1000) integrals 
which can be evaluated with standard numerical meth- 
ods. We have used the Divonne and Cuhre integration 
algorithms from the Cuba library [2^. 

When we apply this method naively to diagrams such 
as the last diagram in Fig. [T] with a sbottom-gluino bub- 
ble inserted in a bottom propagator, the numerical inte- 
gration shows a very poor convergence. Our technique is 
based on a deformation of Feynman parameters: 

Xi ^ Zi — Xi ~i" T'\iyi 

The function xji is constructed using the prescription of 
[^,24], and it depends on Feynman parameters and the 



mass parameters. We can usually set a common numer- 
ical value for the constants Aj, which dictate the magni- 
tude of the deformation, when all mass values are of the 
same order of magnitude. However, the large hierarchy 
of the squared mass values for the bottom quark and the 
gluino or the bottom squarks may yield y^'s of disparate 
magnitude for different Feynman parameters. Small val- 
ues of Ui may not distance the contour of deformation 
sufficiently from singularities. In turn, large yi values 
may result to a contour which is not equivalent with the 
integration over the original Feynman parameters. We 
can remedy these problems by adjusting the parameters 
\i independently for each Feynman parameter. We per- 
form a first Monte-Carlo integration with a small number 
of integrand evaluations setting all to a common value 
A. We use this integration in order to obtain an estimate 
of the maximum values y™"^ of yi. We then compute 
the integrals setting Ai ^ A/y™""^, where A has a typical 
value 0.1 — 0.5. With this refined selection of the param- 
eters which determine the deformation of the integration 
contour we are able to evaluate all diagrams efficiently. 

We were also able to compute the problematic dia- 
grams with a uniform selection of \i values by using a 
Feynman parameterization similar to the one in Ref. [23| . 
This parameterization casts diagrams with bubble sub- 
graphs as a sum of two terms. One of them corresponds 
to an one-loop integral and matches to the counter-term 
for mass renormalization. The second term corresponds 
to a two-loop integral which can be evaluated without a 
specially tuned contour deformation. 

We have used conventional dimensional rcgularization 
(DREG) dimensional reduction (DRED) Jil HI 

and the four dimensional helicity scheme (FDH) [29[ with 
minimal subtraction. In all schemes an anti-commuting 
75 prescription was employed. Similarly to ^1^], we de- 
couple the top-quark, squarks and gluino from the run- 
ning of the strong coupling as. In addition, we perform 
a pole renormalization of all masses and renormalize the 
squark higgs couplings as in [s^]- 

DREG is known to violate the symmetries of the 
MSSM. DRED is a consistent renormalization scheme for 
the computation of the MSSM amplitude through the 
next to leading order in the strong coupling. The two 
schemes are, however, relate d by finite shifts in the cou- 
pling constants and masses [S^]. This equivalence only 
holds when new terms are added to the Lagrangian in 
DRED S^, 33, 13]: operators with e-scalars which arise 
when the 4-dimensional gluon field is split into a D- 
dimensional part and its remaining A — D components. 
Notably, in the MSSM, the only relevant term for the 
amplitude gg h is the mass term for the e-scalars. As 
this mass can be always absorbed into a redefinition of 
the squark masses [sl] , the calculation in the MSSM can 
be performed by setting it to zero. 

The picture described above changes dramatically 
when considering a theory with less symmetries than the 
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MSSM, where supersymmetry is only softly broken and 
the possible e-scalar couplings are still very restricted. 
We observed that, if either SU{2)]^ or (softly broken) 
supersymmetry are absent, a new coupling between the 
e-scalars and the higgs fields emerges radiatively. This 
coupling is indispensable in order to render the one-loop 
cross-section for the process h ^ eeto be of 0{e) (evanes- 
cent). In the MSSM this happens automatically due to 
a cancellation among contributions from up and down 
type quarks and the corresponding squarks. In the SM 
this cancellation does not take place. Most importantly, 
the e-higgs coupling and its consistent renormalization 
should be included in order for the results for the SM 
two-loop amplitude gg ^ h in DRED and DREG to 
agree. 

We compared the results in the two schemes after we 
accounted for known shifts in the strong coupling ag and 
the mass parameters between the two schemes. We found 
that the two results agree only if an additional shift in the 
higgs-squark-squark coupling is performed. We find that 
the relation between the renormalized coupling rn^Vhqiqj 
in DREG and DRED is, 
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where Vkq^q^ is the dimensionless part of the tree 
coupling, and (5^^^^ ,(5^^^*^ are the pole mass- 
renormalization counter-terms in the two schemes. We 
have verified that the shift in Eq. [T] is also needed for 
the computations of the decay rate h — > (jiqj to agree in 
DRED and DREG at one loop. The reason for this shift 
is simple. The one loop corrections to the decay ampli- 
tude h — > Qiqj are identical in DRED and DREG before 
renormalization. However, the renormalization scheme 
for the squark higgs coupling in [30| involves the quark 
pole mass, which is related to the bare mass trough a 
scheme dependent expression. The shift in Eq. [1] simply 
cancels that dependence in such a way that the DRED 
result for the decay rate is recovered. 

We can avoid to compute the contribution of operators 
with e-scalars to the gg ^ h amplitude for the MSSM by 
using the FDH scheme. DRED and FDH treat exter- 
nal polarizations differently; this gives rise to a relative 
factor of (1 — e)^ in the two-schemes for the squared am- 
plitudes through 0{al). Diagrams with internal e scalars 
in DRED are accounted for in the FDH scheme by dia- 
grams with internal Z?s-dimensional gluons, where Ds is 
the dimension of the spin algebra. The dimensionality of 
the loop integrals D is kept distinct with D < Ds; after 
performing the loop-integrations an analytic continuation 
of Ds to 4 dimensions [29| takes place. The FDH scheme 
cannot account for the contributions of diagrams with 
the e-higgs coupling. As we discussed, these are needed 
for the gg ^ h two-loop amplitude in less symmetric the- 



ories than the MSSM, such as the SM. The discrepancy 
cannot be absorbed in any coupling or mass redefinition, 
and the FDH result is inconsistent with the results in 
DRED and DREG for the SM two-loop amphtude. 
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FIG. 2: The squared UV and IR renormalized amplitude for 
a light Higgs boson through Oia'l) and the corresponding K- 
factors 



We now present numerical results for the MSSM two- 
loop amplitudes gg h, H. We neglect here the Higgs 
couplings to quarks and squarks other than the ones in 
the third generation. The two-loop amplitudes are in- 
frared divergent with poles up to second order in the 
dimension parameter e. The singular part is univer- 
sal (soj and cancels against other universal contributions 
at the same order in «<; from real radiation processes. 
We present here the finite part after UV renormalization 
in the MS scheme and subtracting the infrared counter- 
term of Ref. 



37|. 



A complete phenomenological analy- 
sis requires the inclusion of the non-singular parts from 
real radiation and will be the subject of a future pub- 
lication. Nevertheless, useful conclusions may also be 
inferred from solely the two-loop amplitudes since they 
include all diagrams with more than one massive inter- 
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nal particle. These are the diagrams which had not been 
computed earlier in the literature ll|, [iSj, and 



they form a subset which is infrared finite. 

In FigureOwe present our results for the production of 
a light higgs in gluon fusion in the MSSM. On the upper 
panel we show the renormalized squared amplitude for a 
light neutral Higgs boson gg h through order 0{al), 
averaging over external gluon polarizations and colors. 
As discussed above, we subtracted a universal infrared 
counter-term in order to obtain a finite result. We show 
the SM value as a reference. The lower panel shows the 
corresponding K-factor: the ratio of the squared ampli- 
tude through ©(af) divided by the 0{al) result. We also 
include the K-factors obtained in various approximations 
of our full result. We write our results in terms of the 
with a^^{Mz) = 0.1176. We choose a relatively small 
mass value for one of the stop squarks m^^ = 150 GeV 
and vary the mass of the heavy stop quark rrij:^. In the 
bottom sector we have set nih — 5 GeV, rrij^^ — 350 GeV 
and TO^^ = 370 GeV. For the SUSY parameters, we have 
chosen a = 3°, tan/3 = 20, fi — 300 GeV, and the squark 
mixing angles 9^ — — 40°. Finally, mgiuino = 500 GeV, 
mil — 115 GeV, the renormalization scale is fixed to 
the higgs mass, ^^cn = ttih- The renormalization of 
the squark mixing angles is done as in [15] at a scale 
fie = 200 GeV. 

The squared amplitude in the MSSM decreases for 
growing m^^ . This is due to large cancellations between 
diagrams where a top or a light stop couple to the Higgs 
boson. The 0{al) contribution ranges from 15% to 
—40%, becoming negative at large values of m^^ . A signif- 
icant part of this correction originates from the infrared 
finite subset of diagrams with giuinos and squark quartic 
couplings; their contribution to the squared amplitude is 
negative and grows in absolute value with growing mj;^. 
The contribution from bottom and sbottom loops is be- 
low 3% for this value of tan /3, and smaller for lower values 
of tan/3. We also show the K-factor obtained with the ef- 
fective theory calculation of Ref. [l5l| using the published 
program evalcsusy. The approximation given by this 
calculation is remarkably good. At small values of m^^ 
it almost coincides with the MSSM result when neglect- 
ing the bottom contributions whereas at larger values of 
the mass splitting it differs only by a few percent. The 
growing (small) discrepancy is due to the contribution of 
the heavy stop decreasing with a larger mass rrif^ ren- 
dering the contribution of light stop loops, which are not 
perfectly approximated by the effective theory, more sig- 
nificant. 

The scenaria with a light scalar quark have received re- 
cent attention. In particular, a large mass gap in the two 
scalar-t op p hysical states reduces the "fine-tuning" of the 
MSSM [38|; this may be an outcome of supersymmetry 
breaking via "mirage" mediation [s^ . In Fig [2] we ob- 
serve that the squared amplitude decreases as the squark 



mass splitting increases. The parameter region where 
both stop squarks are light is excluded by measurements 
of the p— parameter [3^. However, as the stop mass dif- 
ference increases we obtain viable parameter values; they 
result to a ggh interaction which is significantly weaker 
than in the SM; the discovery of a light Higgs boson may 
then become very difficult at the LHC. 
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FIG. 3: The squared UV and IR renormalized amplitude for 
a heavy neutral Higgs boson through 0{ol'1) and the corre- 
sponding K-factors 



In Figure [3] we show the corresponding results for the 
production of a heavy Higgs. In this case we fixed the 
heavy stop mass m^^ = 350 GeV and varied the heavy 
Higgs boson mass. We set the renormalization scale 
to /^rcn = TiT-H- AH the othcr masses and parameters 
are identical to the light Higgs boson calculation de- 
scribed above. At niH = '2-mi_^ — 300 GeV there is a 
threshold, where the perturbative calculation diverges. 
In the plots, we have conservatively whitewashed a win- 
dow of 5 GeV around the threshold, but we have checked 
that our numerical calculation works fine for phase space 
points much closer to this threshold. We find that the 
©(a^) corrections in the MSSM are very mild, growing to 
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about 20% in the region of the threshold but ammount- 
ing to only a few percent anywhere else. The bottom 
sector contributions are, however, very important, am- 
mounting to almost 40% for small Higgs boson masses. 
As in the light Higgs boson case, contributions from the 
diagrams with gluinos and squark quartic couplings are 
substantial. 

In this paper we have computed the full two-loop am- 
plitudes gg ^ h and gg ^ H m the MSSM, a compli- 
cated extension of the Standard Model. We have devel- 
oped a powerful new method for multi-loop calculations 
and applied it to compute the first two-loop three-point 
amplitude known in the MSSM for arbitrary masses of 
sparticles. Our results will improve the precision of cross- 
section predictions for the gluon-fusion process in the 
MSSM. We are looking forward to further applications 
of our method. 
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